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q \ ABSTRACT 

r/} ' We critically examine how well the evolution of large-scale density perturbations 

d . is followed in cosmological iV-body simulations. We first run a large volume simulation 

and perform a mode-by-mode analysis in three-dimensional Fourier space. We show 
that the growth of large-scale fluctuations significantly deviates from linear theory 
predictions. The deviations are caused by nonlinear coupling with a small number of 
QO ' modes at largest scales owing to finiteness of the simulation volume. We then develop 

\ an analytic model based on second-order perturbation theory to quantify the effect. 

00 ■ Our model accurately reproduces the simulation results. For a single realization, the 

' second-order effect appears typically as "zig-zag" patterns around the linear-theory 

, prediction, which imprints artificial "oscillations" that lie on the real baryon-acoustic 

■ oscillations. Although an ensemble average of a number of realizations approaches the 
00 | linear theory prediction, the dispersions of the realizations remain large even for a 

. large simulation volume of several hundred megaparsecs on a side. For the standard 

ACDM model, the deviations from linear growth rate are as large as 10 percent for a 
• j-h ! simulation volume with L = 500ft." 1 Mpc and for a bin width in wavenumber of Ak = 

■ 0.005/iMpc -1 , which are comparable to the intrinsic variance of Gaussian random 
^ ' realizations. We find that the dispersions scales as oc L~ 3 / 2 Ak~ 1 / 2 and that the mean 

dispersion amplitude can be made smaller than a percent only if we use a very large 
volume of L > 2h~ 1 Gpc. The finite box size effect needs to be appropriately taken 
into account when interpreting results from large-scale structure simulations for future 
dark energy surveys using baryon acoustic oscillations. 

Key words: cosmology:theory - large-scale structure of Universe - methods:N-body 
simulations 
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1 INTRODUCTION 

Understanding the nature of dark energy that dominates the 
energy content of the universe is one of the main challenges 
in cosmology. The time evolution of the mysterious dark 
component is accessible only by astronomical observations. 
Baryon acoustic oscillations (BAO) can be used as a stan- 
dard ruler by which precise meas urement of the cosmological 
dista nce scale is achievable (e.g.. lEisenstein. Hu fc Tegmarkl 
1 19981 ; Seo & Eisenstein 2003; Matsubara 2004). 

Recent large galaxy redshift surveys, the Sloan Dig- 
ital Sky Survey and the 2-degree Field survey, detected 



the signature of the baryon acousti c peaks and thus pro- 
vide constraints on the dark energy (lEisenstein et all 20051 ; 
ICole et al. I l2005l ; iPercivaT et al. I l2007l ; lokumura et al.ll2007f ). 
Future observational programs will utilize the distribution of 
millions of high-redshift galaxies to detect BAO with higher 
accuracy. In order to properly interpret these observations, it 
is necessary to make accurate theoretical predictions for the 
length scale and other characteristic features of BAO (e.; 
iNishimichi et all 120071 ; ISmith. Scoccimarro fc ShetlJ 12008 
Theoretically, a crucial issue is the no n-linear evolution 
of m a tter and galaxy dis t ributions (e.g., Seo fc Eisensteir] 
120051 ; lAngulo et all 120071 ; iGuzik. Bernstein fc Smith! 120071 ; 
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ISmith. Scoccimarro fc Shethl 120071 ). One usually resorts to 
using cosmological TV-body simulations for this, but various 
effects -both physical and numerical- need to be understood 
in order to extract useful information. First of all, the power 
spectrum for a realization of a Gaussian random field has 
intrinsic deviations from expected values at any wavenum- 
ber, i.e., the mode amp litudes are Rayleigh-distributed (see 
e.g., iMatsu bara 2007a|). A realization may thus show an ad- 
ditional oscillatory feature on large scales which compro- 
mises the true BAO signature (Huff et al. 2007). There are 
also a number of numerical issues. Accurate time integration 
is necessary in order to follow the evolution of large-scale 
density perturbations which have small amplitudes. Finite- 
box size limits the sampling of wavenumbers at the largest 
scales, whe re the power amplitud e is dominated by only a 
few modes |Bagla fc Prasadll2006l studied the finite box size 
effect on the mass function of dark matter halos.) 

In this paper, we examine how accurately the evolution 
of large-scale density perturbations is followed in standard 
cosmological TV-body simulations. In particular, we study 
the characteristic "wiggle" features which are often found 
in the matter power spectra calculated from TV-body sim- 
ulations in previous studies. We use an approach based on 
perturbation theory to study nonlinear effects in detail. A 
further extensive study is presented in a separate paper by 
Nishimichi et al. (in preparation). 

Throughout the present paper, we adopt the standard 
ACDM model with matter density Q m = 0.241, baryon den- 
sity fib = 0.041, cosmological constant Ha = 0.759, spectral 
index n s = 0.958, amplitude of fluctuations as = 0.76, and 
expansion rate at the present time Ho = 73.2km s~ Mpc - , 
consistent with the 3-year WMAP results (Spergel et al. 
2007). 
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Figure 1. We plot the evolution of the power spectrum from 
the initial epoch (black line) to z = 3 (green), z = 1 (blue), 
and z = (purple). The measured power spectrum is divided 
by the no-wiggle model of Eisenstein & Hu (1998). We subtract 
the intrinsic deviations from the input power spectrum at the 
initial epoch. The numbers indicate integer sums of n\ + n\ + n| 
of wavenumber vectors. The dashed lines are the one-loop power 
spectra at each redshift (see text). 



2 METHOD 

2.1 The cosmological simulations 

We use the cosmologica l simulation code Gadget-2 
jSpringel. Yoshida fc White] l200ll ; ISpringell [20051 ). For our 
fiducial runs, we employ 256 3 particles in a volume of 
L = 500ft" 1 Mpc on a side. We dump snapshots at a number 
of time steps (redshifts) to study the evolution of the den- 
sity power spectrum. The simulation parameters are chosen 
such that sufficient convergence is achieved in the measured 
power spectrum at the present epoch (Takahashi et al., in 
preparation). 

We generate initial conditions for our runs based on the 
standard Zel'dovich approximation using the matter transfer 
function calculated by C AMB (Code for Anisotropies in th e 
Microwave Background; iLewis. Challinor fc Lasenbvl [2000). 
The initial redshift is set to be Zi n = 30. When we generate 
a realization for a Gaussian random field, the amplitude 
of each fc-mode is assigned such that the ensemble follows 
the Rayleigh distribution. While the mean of the power is 
expected to approach the input value at fc for an ensemble of 
large modes, the actual assigned power in a finite fc-bin can 
deviate significantly from the expected value. Note also that 
a Rayleigh distribution has a positive skew, which causes the 
median to be smaller than the mean. 



2.2 Fourier mode analysis 

We first compute the density field for each output of the 
TV-body simulation. We use the CIC (cloud-in-cell) interpo- 
lation when assigning particles on grids. We check that the 
interpolation method does not affect the scales of interest 
(fc Ss 0.1) by comparing various schemes. We then apply a 
Fast Fourier Transforrr[j to obtain the density field <5(k) in 
three-dimensional Fourier space. We will examine both the 
amplitudes and the phases in detail in subsequent sections. 

In order to study closely the Fourier mode-coupling, we 
calculate the mean amplitude of modes for a given realiza- 
tion with wavenumber vector k = (fci, k2,k%) as 

|k|=* 

where the summation is for all the wavenumbers of |k| = 
fc = (fc 2 + k\ + fc!) 1 / 2 , TVfe is the number of modes in fc, and 
the wavenumber is discretized as ki — (2?r /L)m with an 
integer m. An ensemble average of a number of realizations 
provides its expectation value of P(k) = (P(k)). 

In order to study the evolution of power spectrum, we 
divide the measured power spectrum in equation (fj) at red- 
shift 2 by the initial one at Zi n = 30, and then multiply it by 
the input power spectrum. In this way, the initial random 
scatter included in the power spectrum is removed. 



Simulations of BAO I 



3 



3 RESULTS 

Fig[H shows the evolution of power spectrum P(k) for a 
single realization. We show the mean amplitude for modes 
which have exactly the same wavevector norm, |k| 2 = fcf + 
fcf + fef, rather than binning in k. The vertical axis is the 
power spectrum divided by the no-wiggle model of Eisen- 
stein & Hu (1999). The black line with symbols is the linear 
theory prediction with CAMB. The green, blue, and purple 
lines with dots are the measured mean values at each wave 
number at z = 3, 1, and 0, respectively. The numbers in the 
figure indicate integer sums of n\ + n\ + n\ of wavenumber 
vectors. 

As clearly seen in the figure, the power amplitudes de- 
viate from the linear theory prediction at low redshifts. 
The deviations appear to grow in time monotonically. 
Some modes (e.g. n 2 = 4, 13, 19, 25, 27) grow more rapidly 
than the linear growth, while other modes (e.g. n 2 = 
8, 12, 20, 26, 32) grow less. These features can b e seen even 
in higher resolution simulation of lSpringel at al.1 (|2005l ) (see 
their Fig. 6). Since the initial randomness of the amplitude 
of each mode has been already subtracted in the figure as 
described in section 2.2, the remaining differences plotted 
in Fig. 1 are due either to numerical integration errors or 
to some unknown physical effects. The deviations are in- 
deed large, with the amplitudes being more than 10% at the 
scale of the first-peak of the BAO. It is thus important to 
understand and correct the apparent oscillatory features if 
these are artificial effects. 

In the next section, we show that the deviations are not 
owing to numerical integration errors but due to the finite 
number of modes at the largest scales. We use second-order 
perturbation theory to explain the systematic deviations. 



4 PERTURBATION THEORY 

Second-order perturbation the ory describes the evoluti on of 
a density perturbation as (e.g. iBernardeau et al"1l2002h 



An 



An 



(2) 



where <5i(k) and An are the linear density and the linear 
growth factor evaluated at the initial redshift. The second- 
order term is given by 



S 2 (k) = J2 F * (P. k - PMi (P)*i ( k - P) . 



with 
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We sum up all the modes up to the Nyquist fre- 
quency (256 3 modes in total) in equation (|3}- Here, 
equation (HI includes the fa s test growing mode. 
IBernardeau. Crocce fc Scoccimarrol (|2008h recently present 
the correct formula of A including the sub-leading growing 
mode. 

Let us explicitly write the amplitude and the phase of 
a mode as 

1 FFTW Home page : |http: / /www.fftw.org/| 



*CM) = |*(k,*)|exp(i0(k,z)), (5) 
Then the evolution of amplitude in each mode is 
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whereas the phase evolution is 
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up to second order. The expressions in equations (J6]l and 
(0 are independent of the initial redshift for the late time 
(D ^> An), since 81 oc An and 82 oc P(k, Zi n ) oc An- We 
do not distinguish between 8 and Si at the initial redshift 
(zin = 30), since 82 is much smaller than <5i at that time. Q 
and provide more detail analysis. 

It is clear from equation @ that nonlinear mode- 
coupling occurs with particular sets of wavenumber vectors 
such that p + q = k. From equation @, we obtain 



A(p,k-p) 

for k <C p, and 
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A(p,k-p) -> --cosf 
2 p 

for k ^> p. Here 9 is an angle between k and p. Hence 
the coupling to the mode of much smaller scale p(^> k) is 
negligibly weak, while the coupling to much larger scale p(<^i 
k) is strong. In summary, most of the contribution to the 
second-order evolution of a mode comes from the modes of 
comparable scales or larger[f| 

For a Gaussian random field, the mode amplitudes are 
Rayleigh-distributed, and thus there is a finite probability 
that a mode has a very large or a very small amplitude with 
respect to the expected mean value. Some peculiar modes, 
which have very large or very small amplitudes compared to 
the mean, strongly affect the growth of other modes through 
the mode-coupling as described in the above. 

In an ideal situation where there are infinite number of 
modes, the second term in equation © vanishes. In that 
case, the leading correction arises from the forth order of 
Si. Then the resultant power spectrum with the one-loop 
correction is, 



Aloop(fc, z) = 



An 



Pu(k)+ 



DM 

An 



[P 22 (fc)+Pi 3 (fc)],(10) 



where Pu = (\Si\ 2 ). P 22 = {lfe| 2 V Pi 3 = KRefeja 
ilMakino. Sasaki fc Sutol Il992l ; | Jain fc Bertschingerl Il99 
iJeong fc Komatsul 12009 ). We integrate from k — 2n/L to 
the Nyquist frequency in the calculation of P22 and P13. 



2 Nishimichi et al. (in preparation) distinguish S from 
i 5i at the initial epoch with the 2LPT initial condition 
llCrocce. Pueblas fc Scoccimarroll2006h 

3 iMuecket et al.lll98sT examined the growth of the small-scale 
perturbation on the background of the large-scale perturbation. 
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Figure 2. Evolution of the deviation of the power amplitude with respect to the linear theory prediction. The dots are the measurements 
from our simulation, and red solid lines are the model prediction using the second-order perturbation theory. The integers denote 
n 2 = n 2 + n?, + n| of wavenumbers, and the figures show different range of n 2 , n 2 = 1 — 8 (upper left panel), n 2 = 9 — 16 (upper right 
panel), n 2 = 17 — 24 (lower left panel), and n 2 = 25 — 32 (lower right panel). 



The dashed lines in Fig[T]are the one-loop power spec- 
trum at each redshift. It suggests that the linear theory is 
applicable for k < 0.07/i/Mpc at z — 0. However, the finite 
mode coupling in the second term of equation signifi- 
cantly changes the evolution of the power spectrum even in 
the linear regime^ 

FigOS shows the evolution of the mean amplitude of 
modes with identical wavenumber n 2 in the range of 1 — 32. 
Here, n 2 ~ 30 corresponds to the position of the first peak 
(see Fig[T}. The four panels are for n 2 = 1 — 8 (upper left 
panel), n 2 = 9 — 16 (upper right panel), n 2 = 17— 24 (lower 
left panel), and n 2 = 25 — 32 (lower right panel). The dots 
are the measurement from simulation outputs, and red solid 
lines are the theoretical prediction from the initial density 
fields at Zi n — 30 in equation (B). The second-order pertur- 
bation theory reproduces the simulation results rather well. 

4 ISetd dl999T> also investigated the finite mode effect on the one- 
loop correction terms, P22 + P13, in equation (I T D . 



The theory fits the data within 0.5% at z = 2 and 2% at 
z = for larger scale (n 2 = 1 — 8), whereas within 1% at 
z = 2 and 10% at z = for smaller scale (n 2 = 25-32). This 
is because the second order perturbation theory is applicable 
at large scales and/or at high redshift. 

Fig'EJis the same as FigO but for phase evolution. We 
plot the results only for modes with m ^ 712 ^ H3, because 
the mean of the phase at k, ^ 0(k), is zero (since (/>(k) + 
cj>{— k) = 0). The phase shifts are typically w 0.1 radian at 
z = 0. Perturbation theory well reproduces the results. Even 
if there are infinite modes, the right hand side of equation 
J7J still remains. The phase shift is not due to the finite box 
size effect. 

Previously iRvden fc Gramannl (|l99lf ) and iGramarml 

( 1992) studied the evolution of amplitude and phase in each 
mode using two dimensional simulations. They also calcu- 
lated second-order perturbation theory and found the de- 
viation from the linear theory grows in proportiona l to th e 
scale factor in the EdS model. ISuginohara fc Sutol l|l99ll) . 



Simulations of BAO I 5 




ISoda fc Sutol (|1992) and lJain fc Bertschinger] (|l998l ) also ex- 
amined the nonlinear evolution in each mode. However they 
did not compare the theoretical prediction with the simula- 
tion results in detail. Their motivations were to understand 
the evolution of the density fluctuations in the nonlinear 
regime, whereas our interest here is in the growth of pertur- 
bations at the linear scale. 



5 STATISTICAL ANALYSIS 

The previous section considers second-order effects for a sin- 
gle realization. In this section we run 100 simulations to cal- 
culate dispersions of amplitude and phase deviations from 
linear theory. We prepare the 100 realizations for each of 
three box sizes of L = 500/i _1 Mpc, lfo _1 Gpc, and 2/i _1 Gpc, 
and Zin = 30, 20 and 10, respectively. 

Fig(4] shows the remaining amplitude dispersions from 
the linear theory prediction after correcting for the initial 
randomness at z — for L = 500/i _1 Mpc (top), l/i -1 Gpc 
(middle), and 2/i _1 Gpc (bottom). Since we already sub- 
tract the initial deviations due to the Gaussian distribution, 



the residuals arise from the mode-coupling during the evo- 
lution. The grey dots with error bars are the means with 
la scatters. By using a sufficiently large number of realiza- 
tions, the means converge to the true values (solid line), and 
the magnitude of the dispersions is insensitive to the num- 
ber of realizations. For L — 500/i _1 Mpc, the dispersions are 
~ 10% near the first peak, and ~ 5% even for a very large 
volume of 2/i _1 Gpc on a side. The dashed lines show the 
theoretical prediction of the la scatter, which is the rms 
(root-mean-square) of the second term in equation (|6]) : 

/ / p(k,z)/p(k, Zin ) y\ 

amp " \V D{ z y/Di l ) I 

4P 22 (fc, Zin ) 1 (D{z)\ 2 

Pll(Min) AiV fc \D m J ' 1 ' 

Here, AA^ is the number of modes in the bin, ANk = 
Aim 2 An with n = (L/2w)k. In this unbinning case, the num- 
ber of modes is AiVfc = kLAn 2 (with An 2 — 1). The dashed 
lines well reproduce the results. 

Fig(4]also shows the results for the binned data of Ak = 
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Figure 4. The amplitude dispersions of the 100 realizations at 
z = for L = 500/1" 1 Mpc (top), lh" 1 Gpc (middle), and 2ft- 1 
Gpc (bottom). The grey dots with error bars are for the un- 
binned data, while the black big symbols are for the binned data 
of Afc = 0.005/i/Mpc. The value of k for the binned data is the 
weighted mean of k with the number of wavenumbers in the bin. 
The dashed lines are the theoretical prediction. 
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Figure 5. The amplitude dispersions calculated from our simula- 
tion outputs (filled circle •) and the theoretical predictions (solid 
lines). We also show the dispersions due to the initial Gaussian 
distribution (dashed lines). The vertical dotted line is the position 
of the BAO first peak. 
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Figure 6. We compare two dispersions. Blue points with error 
bars show intrinsic scatter around the expected mean power spec- 
trum for initial Gaussian random density fields. Black points show 
the dispersions owing to the finite nonlinear mode-coupling effect. 



0.005/i/Mpc by the black big symbols. In this case, we use 
the power spectrum defined as P(k) = (1/ANk) |<5(k) | 2 , 
summing up all the modes between (k — Afc/2, k + Afc/2), 
instead of equation {l}. Here, the number of modes in the 
bin is 

AN k = (L 3 fc 2 )/(27r 2 )Afc. (12) 

We calculate the means and error bars for the binned P(k). 

Fig[S] shows the amplitude dispersions calculated from 
our simulation outputs for Afc = 0.005/i/Mpc (filled circle) 
and the theoretical prediction (solid line). From this figure 
with equations (|lip and (|12p . we find that the dispersion is 
approximated as 

{Z = 0) ^ 2% (iG^Th) (cIoC^M " 1/2 ' (13) 
at k = 0.02 — O.lft/Mpc. The dispersion is proportional to 
AN~ 1/2 oc L~ 3/2 Ak~ 1/2 from equation fl2]l. Note that even 



with a large simulation volume of L ~ 1 Gpc with fc-binning, 
the dispersions still remain at the level of a few percent. 

So far we have discussed the amplitude of deviations 
from linear theory. Here we also consider the intrinsic scat- 
ter of the initial Gaussian random realizations. In Fig[5] 
the dashed line is the dispersion for the initial distribu- 
tion, which is given byQ (AN k /2y 1/2 . Fig[S]shows that the 
dashed lines decrease as oc (AAT fe )- 1/2 oc AT 1 , while the solid 
lines increase because P22/P11 increases (see equation |11]). 
These two dispersions are comparable at k ~ 0.1ft/Mpc 
where 2P22/-P11 — 1 at z = 0. About a half of the dispersions 
near the position of the BAO first peak (k ~ 0.07/i/Mpc) 
are attributed to the second-order effects. The result sug- 
gests that, at large scales, k < 0.1/i/Mpc, the dispersions 
arise mainly from the initial Gaussian random distribution, 

5 The number of modes ATVj. is divided by 2 because the Fourier 
modes of 5(k) and S(— k) are not independent. 
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Figure 7. The phase dispersion of the 100 realizations. The solid 
line is the theoretical prediction. 

while at smaller scale k > 0.1/i/Mpc they are from the mode- 
coupling (based on the second or higher order perturbation) 
during the evolution. In Fig(6]the blue symbols are the re- 
sults for our 100 realizations. The black symbols are same as 
in the top panel of Fig|4]for Afc = 0.005/i/Mpc. As expected, 
the initial random realizations have larger scatters around 
the mean expected power spectrum, especially at the largest 
scales. 

We have also performed a similar analysis for the evolu- 
tion of the mode phases (equation [7]). Fig0shows the phase 
dispersion calculated from our simulations (the dots). Here 
we set — tv ^ (<fi — <f>i n ) ^ 7r and calculate (|<5i| 4 (</> — </!>in) 2 } 
instead of {(cj> — cj>i n ) 2 }. This is because (<f> — (j>i) oc 1/Si in 
Eq.(0 and its dispersion diverges at Si = 0. We obtain the 
phase dispersion from equation (0 as[f| 

(|Mk)| 4 [0(M)-<Mk)] 2 ) = P 22 (fc,*in) [D(Z)Y 

(|<5i(k)| 4 ) 6Pn(fc,z in ) ^ An ) A ' 

The solid lines are the theoretical prediction, which fit the 
simulation results well. The phase dispersion in equation 
<|14[) . as well as the amplitude dispersion in equation (111[) . are 
independent of the initial redshift. In the non- linear limit of 
k — > oo, the phases are distributed rand omly, and the phase 
dispe rsion approaches to tt/^/S rad (e.g. lRvden fc Gramannl 
S). 



6 DISCUSSION AND CONCLUSIONS 

In this paper, we critically examined how accurately cosmo- 
logical iV-body simulations describe the evolution of large- 
scale density distributions, particularly focusing on the lin- 
ear and/or quasi-linear scales. For the power spectrum cal- 
culated from a single realization, we found that the growth of 
large-scale fluctuations significantly deviates from the linear 
theory prediction, and the enhanced or suppressed growth 

6 Ijain &: Be rtschingerl lll996f ) previously derived equation 1141 1 
with an approximation for the long-wave mode coupling. 
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Figure 8. The power spectrum at z = 2. The green line is the 
simulation output. In the red line, we subtract the second-order 
perturbation contribution from the simulation output. The blue 
line is the one-loop power spectrum. 



of perturbations produces an ugly noisy pattern in the mat- 
ter power spectrum. This deviation is not due to the nu- 
merical errors in the TV-body code, but due to the non- 
linear coupling between finite numbers of modes originat- 
ing from the finite box size. To study the effect of the fi- 
nite mode-coupling in detail, we developed perturbation the- 
ory and quantitatively estimated the finite-mode coupling to 
the power spectrum amplitude. Mode-by-mode analysis in 
three-dimensional Fourier space reveals that the finite mode- 
coupling from the second-order perturbation is sufficient to 
explain the deviation from linear theory prediction on large 
scales. The dispersion of the mode-coupling effects estimated 
from second-order perturbation scales as oc L~ 3 ^ 2 Afc -1 ^ 2 , 
and this may surpass the intrinsic scatter of the initial Gaus- 
sian distribution. Since the finite mode-coupling does not 
vanish even for a large-volume simulation, it is of critical 
importance to correct it properly for high-precision studies 
of baryon acoustic oscillations. 

We show that the perturbative approach is very helpful 
to quantify the significance of finite-mode coupling and this 
can be utilized as an efficient and powerful tool to correct the 
finite-mode coupling. As an example, in Fig. [8] we evaluate 
the power spectrum directly obtained from a single realiza- 
tion at z = 2, and subtract the finite-mode coupling using 
the second-order perturbation. Compared the result before 
subtraction with that after subtraction, the deviation from 
linear theory is dramatically reduced and the noisy struc- 
tures are effectively wiped out. As a result, even the single 
realization data of iV-body simulation faithfully reproduces 
the linear theory prediction on large scales. 

Although the present paper mainly concerns with the 
second-order perturbation theory, higher-order perturba- 
tions are also important for the relevant scales of the mea- 
surement of baryon acoustic oscillations, where the acoustic 
signature te nds to be erased by the effect of non-linear clus- 
tering (e.g. iMcDonaldl 120071 ; ICrocce fc Scoccimarrol 120071 ; 
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iMatsubaral l2007bl ; iTaruva fc HiramatsvJl2008h . The height 
of the first peak is found to be reduced about 2% (J. Wang, 
A. Szalay et al. in preparation). Thus, the inclusion of the 
higher-order terms may be important for the estimation of 
the finite-mode coupling, which would be helpful to further 
reduce the noisy structures on small scales. 

We note that the variance of the growth of matter 
power spectrum with respect to the linear theory pre- 
diction, {[(P/P in )/(D/DiJ) 2 - l] 2 }, which we have stud- 
ied, is different from the variance of the power spectrum 
itself, ((P - P) 2 ). It remains unclear if the numerical 
effects studied here are important in evaluating covari- 
ance matrices Ce.g.. pcoccimarro. Zaldarriaga fc Hull 19991 ; 



iMeiksin fc Whitdfl999t iNevrinck fc Szapudill2007f ) In future 
work, we will study nonlinear and numerical effects in the 
power spectrum covariance using a large set of simulations 
and analytic models. 
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